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Abstract 

The turbulence module recently developed for 
the NPARC code has been extended to include wall 
functions. The Van Driest transformation is used so 
that the wall functions can be applied to both in- 
compressible and compressible flows. The module 
is equipped with three two-equation K-e turbulence 
models: Chien, Shih-Lumley and CMOTT models. 
Details of the wall functions as well as their numer- 
ical implementation are reported. It is shown that 
the inappropriate artificial viscosity in the near-wall 
region has a big influence on the solution of the wall 
function approach. A simple way to eliminate this 
influence is proposed, which gives satisfactory re- 
sults during the code validation. The module can 
be easily linked to the NPARC code for practical 
applications. 

1. Introduction 

A turbulence module (Zhu and Shih, 1995) has 
been developed for the NPARC code. The mod- 
ule is written in a self-contained manner so that the 
user can use any turbulence model built in the mod- 
ule without concern as to how it is implemented 
and solved. The input to the module is the mean 
flow variables, boundary and geometric information 
which are to be provided by the NPARC code. The 
output of the module is the turbulent eddy- viscosity 
/i t and/or relevant turbulent source terms which are 
needed for the mean flow calculation. The interac- 
tion between the NPARC code and the turbulence 
module will give the final turbulent flow solution. 
The turbulence module is conceived as a vehicle to 
facilitate the technology transfer from the model de- 
velopment to the model application. With the mod- 
ule, any development in the turbulence modeling can 
be quickly and directly validated by and applied to 
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calculations of flows of engineering interests. 

So far, three low Reynolds number two-equation 
turbulence models have been built into the module: 
Chien (1982), Shih and Lumley (1993) and CMOTT 
realisable models (Yang et al., 1995; Zhu and Shih, 
1995). Unlike the Chien’s model, both the Shih- 
Lumley and CMOTT models do not involve the di- 
mensionless wall distance y + , an advantage for sep- 
arated flow calculations. The module with these 
turbulence models has been applied to a number of 
flows including flows over a flat plate, in an ejector 
nozzle, in a transonic diffuser, and a boat-tail noz- 
zle flow (Yang et al., 1995). For all the flow cases 
tested so far, it has been found that both the Shih- 
Lumley and CMOTT models produce improved or 
similar predictions compared with the Chien model. 
The CMOTT model with variable C M turns out to 
be more computationally robust than the other two. 
It was able to give numerical solutions in cases where 
the models with constant C h suffered from numeri- 
cal instability. 

The major problems or difficulties associated 
with the low Reynolds number turbulence models 
are: 1) They require very fine grid spacing in near- 
wall regions, thus increasing considerably computa- 
tional burden, especially in three-dimensional cases. 
Moreover, the highly stretched nature of mesh dis- 
tribution may have an adverse impact on numerical 
stability. 2) Most of models are not of tensorial in- 
variant form, that is, they contain a distance param- 
eter normal to the wall. The wall-distance depen- 
dency causes inconvenience for model applications 
in complicated geometries. Currently, great effort 
is being given in the area of near- wall turbulence 
modeling to remove this dependency, but no satis- 
factory result has been obtained yet. 3) Most of the 
low Reynolds number turbulence models were fine- 
tuned against attached flows, which is, of course, 
not sufficient to guarantee their good performance 
for separated flows. 

An alternative is to use the high Reynolds num- 
ber turbulence models. Here, the governing equa- 
tions are integrated to a point far outside the vis- 
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cons sublayer rather than down to the wall, and 
the near-wall region is bridged over with the wall 
functions. Although in theory the wall function ap- 
proach is only valid for certain attached flows with 
no pressure gradient and mass transfer, it has been 
applied in practice to many separated flows with var- 
ing degree of success. For those flows where maxi- 
mum shear stresses occur far away from the wall, 
the near-wall turbulence modeling is not critical for 
overall flow simulations. In these cases, use of the 
wall functions has a very beneficial effect on the sta- 
bility and economy of computations. Although the 
principle argument for originally adopting the wall 
function approach (economy of grid points) has been 
weakened as larger and faster computers have be- 
come available, it will still find its applications in 
predicting complex flows, especially for large scale 
engineering problems. 

In the present work, we extend the turbulence 
module by including the wall functions. For incom- 
pressible flows, the universal law of the wall may be 
expressed as 


rates Equation (1) into the Navier Stokes equations. 
In this way, there is no need to solve Equation (1) for 
U T by sub-iteration. The implicit method turns out 
to be more stable than the explicit one. Another im- 
portant issue is the artificial viscosity. Chitsomboon 
(1994) found that the artificial viscosity originally 
implemented in the NPARC code totally spoiled the 
solution of the wall functions. This was because the 
artificial viscosity became unrealistically large in the 
vicinity of walls due to very steep velocity gradients 
resulting from the coarseness of grid spacing as re- 
quired by the wall function approach. He fixed up 
this problem by extrapolating velocities at the wall 
rather than using the physical values of no-slip ve- 
locities, when calculating the artificial viscosity. In 
the present work, we simply turn off the artificial 
viscosity in the near- wall region. 

In what follows, we will present the details of 
the wall functions and their implementation, and 
demonstrate how to link the turbulence module to 
the NPARC code. Detailed applications will be re- 
ported in another paper. 



where k = 0.41 and C = 5.2. The derivation of 
Equation (1) is based on the assumption that the 
shear stress in the region close to the wall is constant 
and equal to the wall shear stress. It has been shown 
(Viegas et al., 1985; Huang and Coakley, 1993) that 
the same form also exists for compressible flows with 
the velocity U being replaced by the Van Driest 
transformed velocity U c (Van Driest, 1951). For the 
K-e turbulence models, the convection and diffusion 
terms of their transport equations are negligible in 
the inertial sublayer so that local equilibrium pre- 
vails, which implies that the production of the tur- 
bulent kinetic energy K is equal to the dissipation 
rate e of K. The local equilibrium condition leads to 
two simple relations which can be used as boundary 
conditions for K and e for both incompressible and 
compressible flows. The compressible wall functions 
have been successfully applied to both attached and 
separated flows under Mach number ranging from 
0.1 to 10 (Huang and Coakley, 1993). 

Although the wall functions look simple, their 
numerical implementation is not trivial. The main 
difficulty comes from the logarithmic law in which 
both U and JJ T are unknown, and U T cannot ex- 
plicitly be solved for. It is prone to being numer- 
ically unstable if one uses Equation (1) and itera- 
tively solves U r to obtain the boundary conditions 
for the Navier-Stokes equations. In this work, we 
use an implicit procedure which directly incorpo- 


2. Calculation Approach 

The following presentation is only restricted to 
the numerical aspects related to the implementation 
of the wall functions. Refer to Cooper and Sirbaugh 
(1989 and 1990) about the details of the NPARC 
code and to Zhu and Shih (1995) about the details 
of the turbulence module. 

2.1 Governing Equations 

In the NPARC code, the following nondimen- 
sional, Reynolds-averaged Navier-Stokes equations 
are solved: 


dQ dF x dF 2 dG x dG 2 

dt + d$ + drj ~ ~df + ~d^ ^ ^ 

where Q, Fj and Gj (j = 1, 2) are the conservation 
variable vector, inviscid flux vectors and viscous flux 
vectors, respectively. Because only the viscous flux 
vectors need to be modified with the use of the wall 
functions, their detailed forms are given below: 
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In the turbulence module developed for the 
NPARC code, the three K-e turbulence models 
have been implemented: Chien, Shih-Lumley and 
CMOTT models. The detailed forms of these mod- 
els are given in Zhu and Shih (1995). 

2.2 Wall Functions 

The compressible law of the wall (Huang and 
Coakley, 1993) is used in the turbulence module. 
Following the NPARC nondimensionalization, this 
law can be written as 


v} = W r =\ HEy±] 


( 10 ) 


where 


U T = y/(rjp)71ii 
y + - R t U T y{p/ y.) w a u 
k = 0.41, E = 8.4317 


(11) 


and U c is the Van Driest transformed velocity de- 
fined by (Van Driest, 1951): 


u c = Vb 


sin 


where 


■■ft 2 )-*-®] 
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2T wa ii 
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(13) 


D = VA* + B 


In the near- wall region, with the convection ne- 
glected the energy equation can be reduced to give 
an expression for the total heat flux 


q = qwaii + Ut 


(14) 


and with the local equilibrium assumption (Launder 
and Spalding, 1974), the turbulent kinetic energy K 
and its dissipation rate e can be calculated by 


K = 


Twall/ P 


0.3 


c _ (Wp)*/ 2 

*y 


(15) 


(16) 


In the above expressions, the subscript wall 
refers to the value of the corresponding function at 
the wall. Equations (10) and (14) - (16) form the 
wall functions which are used to bridge over the first 
grid point and the solid wall. 

2.3 Numerical Implementation 

From Equations (10) and (11), the near- wall 
shear stress can be written as 


U c 


T = T wall — PwallUr—r = — (17) 


— — 

’ U? R e y 

where pt is an effective turbulent viscosity connect- 
ing the wall and the first grid point 


p t = 


V+PwallUc 

uu+ 


(18) 


An advantage of Equation (17) in calculating sepa- 
rated flows is worthy of note: the direction of the 
wall shear stress T wa zz is determined by that of the 
flow velocity U while r V9 u calculated from Equa- 
tion (10) or (1) can only have a positive sign. 

From Equation (17), the general form of shear 
stress at the first grid point can be expressed as 


T — T w a n — A Ut 


(19) 


where 
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From geometrical consideration, we have 


f iit/(R e An) if y+ > 11.6 
L A *w<xJz/(£eAn) otherwise 


( 20 ) 


nl + + r%) _ A I 2 _ Ai 

J J-i ~ AO “ An 


17* and An are the tangential component of the re- 
sultant velocity and the normal distance from the 
wall, respectively. Equations (19) and (20) sim- 
ply treat the near-wall region as a laminar sublayer 
(y+ < 11.6) and a fully turbulent layer (y+ > 11.6). 
This treatment prevents the wall function procedure 
from producing abnormal results when y+ tends to 
zero, such as in the vicinity of separation or reat- 
tachment points. 

Similarly, the total heat flux at the first grid 
point can be written as 

q = -a(T-T vall ) (21) 

where 

r A/[( T -1 )P rt ] if y + > 11.6 
a = { (22) 

l A/[(7 - 1 )P T ] otherwise 


where AQ and A / are the volume and face area of 
the control volume, respectively. 

In the wall function approach, all the stresses 
acting on the cell face considered are replaced by the 
wall shear stress given by Equation (19). Therefore, 
for the south wall, Equations (25) are replaced by 


922,$ = AlXU tx 

<723,# = AIXUty 

924,$ = U 922,$ + V <723,* 

+A la(T-T w ) 


(27) 


where U\ x and Ut y are the x- and y-components of 
the tangential velocity 17*, respectively. If n x and 
Tiy are the Cartesian components of the unit normal 
vector at the wall, U tz and Ut y can be calculated by 


and the heat flux at the wall can be calculated by 


Q -wall — — Ol(T Tyjall) Ut * Tuan (23) 

Consider first the wall of 77=constant. In this 
case, only the viscous flux vector G 2 in Equation (2) 
needs to be modified with the use of the wall func- 
tions. In the NPARC code, it is calculated by 


dg22 

drj 


= <722, n ” <722,* 


£<723 

8t) 


= 92Z,n — <723,# 


(24) 


£<724 

= g24,n “ <724,# 

or) 

where the subscripts s and n refer to the south and 
north faces of the control volume under considera- 
tion (Figure 1). From Equations (5) - (7), the com- 
ponents of the vector G 2 can be written as 


<722 


<723 
9 24 


Mtot &£ 

J R t dr) + ' ' ' 

vZ + fftot dv 
J R t dr)*'" 


Ut 21 + Vt 2 2 

+ vlp_JT + 

J 07) 


(25) 


U tx = n\V - n x nyV 
Ut y = n*V - n x iiyU 
Similarly, for the north wall, we have 


(28) 


922, n = -AlXU tz 
92Z,n — -AlXU ty 

924,n = ^<722,* + ^<723,# 

-A l*(T-T w ) 


(29) 


For the wall of £=constant, only the viscous flux 
vector Gi in Equation (2) needs to be modified when 
using the wall functions. In the NPARC code, we 
have 


£<712 

= 9l2,e — 9 12, w 


— gi3,c — 9iz,w (30) 


£<714 

— gl4,« gl4,ti? 

where the subscripts w and e refer to the west and 
east faces of the control volume under consideration 
(Figure 1). From Equations (4), (6) and (7), the 
components of the vector G\ can be written as 
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013 


H + HtHctdU 
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il+zh itotdv 

J R t dt 


9u = Urn +V f 12 
+ J d£ 

and the geometric consideration leads to 


(31) 


g + $ __ J~*(g + ( 2 ) _ Al 2 _ A l 

J ~ jr-» An An { ’ 

where AQ and Al are the volume and face area of 
the control volume, respectively. 

After replacing all the stresses acting on the cell 
face considered by the wall shear stress given by 
Equation (19), Equations (31) become: 
for the west wall 


012,w 

= AlXU tx 

01S,to 

= AlXUty 

014, w 

— U 012, to + V 013, to 


+A la(T - Tuan) 

for the east wall 


012, e 

= -AlXU tx 

01S,e 

= -AlXUty 

014, e 

= 012,0 + ^013, e 


Al&[T 2«rali) 


(33) 


(34) 


The sequence in which the above equations are 
solved together with the Navier-Stokes and turbu- 
lence equations in the code is as follows: 


a. Initialize all field values. 


b. Calculate r wa n and y + using Equations (19) 

and (11). 

c. Fix the values of itT and e at the first grid points 

using Equations (15) and (16). Solve the tur- 
bulence equations. 

d. Calculate q wa u using Equation (23). 

e. Calculate U c and U+ using Equations (12) 

and (10). 

f. Update using Equation (18). 


g. Update a using Equation (22). 

h. Update < 722 , 02 s 024 using Equation (27) 

or (29); or update gi 2 , 013 and g 14 using Equa- 
tion (33) or (34). Solve the Navier-Stokes equa- 
tions. 

i. Return to step b with updated field values. 

The sequence of steps b to i are repeated until the 
calculation converges. 

Note that by definition, the turbulent eddy- 
viscosity fit is zero at the wall, such as in the case of 
the low Reynolds number turbulence models. When 
using the wall functions, Equation (18) introduces 
the effective turbulent viscosity which is defined at 
the wall in the turbulence module. Therefore in 
post-processing, the wall friction coefficient Cj can 
be calculated in the same way as for laminar flows, 
except that the molecular viscosity \i is replaced by 
the turbulent viscosity fi t at the wall. 

3. Module Usage 

The present turbulence module (version 2.0) is 
written based on the NPARC version 2.2. The fol- 
lowing are those parts of the code which may require 
user’s attention when using it. 

3.1 Module 

To facilitate identification, all the subroutine 
names in the module start with CM. In order to 
use the module, the user only needs to call its three 
subroutines: CMAO, CMALL and CMVRHS, in the 
NPARC code. 

Subroutine CMAO. This subroutine transfers 
from NPARC to the module the parameters to define 
flow, geometric and boundary conditions. In addi- 
tion, it has the following user-specified parameters: 

FDEFER — Blending factor in the convection 
scheme. Its value may vary from 0 to 1 with the 
limiting value 0 for the first-order accurate upwind 
and 1 for the second-order accurate HLPA scheme. 
The solution tends to be more stable, but also more 
diffusive when this factor is reduced. 

JTDMA, KTDMA — The value of these integer 
parameters determines whether to use the TDM A 
solution algorithm in J- or K-direction (1 does; 0 
doesn’t). Currently, they are set to 
JTDMA=1 
KTDMA=1 

BDMAX(i), BDMIN(i) — Upper and lower bounds 
for the values of K (i=l), e (i=2) and lit (i=3). 
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These bounds are introduced for numerical purposes 
only, that is, to prevent the corresponding turbu- 
lence quantities from becoming negative or abnor- 
mally large during the solution process. Currently, 
they are set to 

BDMAX(l)=1.0E+6 

BDMAX(2)=1.0E+6 

BDMAX(3)=1.0E+4 

BDMIN(l)=1.0E-8 

BDMIN(2)=1.0E-8 

BDMIN(3)=1.0E-3 

which should cover a wide range of the physically 
meaningful values of K , e and ft*. It is to be 
noted that these values are only valid for the non- 
dimensional turbulence quantities, as defined in the 
NPARC code. 

RELAX(i) — Under-relaxation factors for K (i=l) 
and € (i=2). Currently, they are set to 
RELAX(1)=0.8 
RELAX(2)=0.8 

Should instability occur, try to reduce these values. 

Subroutine CMALL. This is the main subrou- 
tine to control the solution sequence in the module. 
The array variable VIST is the turbulent viscosity 
which is needed in NPARC for calculating tur- 
bulent flows. The array variables TE, ED and YPS 
are K y € and y+, respectively, which can be used for 
post-processing. Normally, there is no need for user 
to change this subroutine. 

Subroutine CMVRHS. This subroutine which 
is the counterpart of the subroutine VISRHS in 
NPARC is for introducing the wall function modi- 
fications into the right-hand side viscous flux terms. 
There is no need for user to change this subroutine. 

3.2 NPARC 

The authors have made all the modifications 
necessary for NPARC to use the module. The follow- 
ing shows where these modifications are in NPARC. 
All the alterations are marked between C<< and 
C>> in the code. 

Namelist TURBIN. The integer parameter 
IMUTR2 is used to select the turbulence models in 
the module with 

IMUTR2=1G1 Chien model 

102 Shih-Lumley model 

103 CMOTT model 

A new integer parameter MWALF is introduced to 
select the near-wall approach with 


MWALF =0 low Reynolds number approach 
1 wall function approach 

Correspondingly, a new statement is added 
in the include file NPARC.INC: 

COMMON /CMOTT/MWALF 
and in the subroutine TURBIN: 

CALL NLGETI(’MWALF\ MWALF). 

Subroutines FILT1, FILT2, FILTER. An array 
FAV01 has been introduced into each of these sub- 
routines to eliminate the artificial viscosity in the 
near- wall region. 

Subroutines INITIA, WREST. The model 
identifier ITURB for each K-c turbulence model in 
the module is given an integer value greater than 
100. To reflect this expanded choice for turbulence 
models, the read and write statements for the turbu- 
lent quantities in these two subroutines are modified 
as 

IF(ITURB.EQ.4 .OR. 

& ITURB.EQ.5 .OR. 

& ITURB.EQ.7 .OR. 

& ITURB.GT.100) 

& READ(2) or WRITE(4) (values of K, e, /**) 

Subroutine MUTURB. The subroutines CM A0 
and CMALL of the module are called here. 

Subroutine STPF2D. The subroutine CMVRHS 
of the module is called here. 

4. Application 

Turbulent boundary layer flow over a flat plate 
with zero pressure gradient was selected as the first 
test case for code validation. The solution of the wall 
function approach (WF) was compared with both 
the experimental data (Exp) of Wieghardt and Till- 
mann (1951) and the solution of the low Reynolds 
number approach (LR). The Shih-Lumley model was 
used in this test. Note that the high Reynolds num- 
ber version of this model is of the same form as the 
standard K-e model of Launder and Spalding (1974). 
Figure 2 shows the flow geometry and boundary con- 
ditions used in the calculation. For the wall func- 
tion approach, grid points were 111x55, the first 
grid points above the wall had the y+ value of 60 
and 14 grid points in the x-direction were located 
before the leading edge of the flat plate. For the 
low Reynolds number approach, grid points were 
111x81 with the distribution of x-points being the 
same as in the wall function case and the first y + 
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being 0.3. Since the NPARC code is for compress- 
ible flows while the experiment ( Wieghardt and Till- 
mann, 1951) to be compared was for incompress- 
ible flows, a freestream Mach number of 0.2 was 
chosen. The reference Reynolds number R t was 
6.9 x 10 6 . All calculations started from an initial 
field with U = 0.2, V = 0,/x t = 1, jBT = 0.005 and 
e = 0.09R e pK 2 /ii t . The influence of the artificial 
viscosity originally implemented in the NPARC code 
was first examined. This was done by simply reduc- 
ing the coefficient DIS4 of the 4th-order artificial 
viscosity. For this particular flow, the 2nd-order ar- 
tificial viscosity was inactive due to the zero pres- 
sure gradient condition. The results are shown in 
Figure 3 for a streamwise mean velocity profile at 
a given location and in Figure 4 for the wall fric- 
tion coefficient. It can be seen that the artificial 
viscosity with DIS4=0.6 totally spoiled the solution 
of the wall function approach and very good results 
were obtained only when DIS4 was reduced to 0.001. 
Although reducing the artificial viscosity can signif- 
icantly improve the prediction in this test case, this 
method is generally unfeasible, because the artificial 
viscosity will also overly be reduced in places where 
certain amount is needed. In the near- wall region, 
the turbulent viscosity prevails, which is enough to 
stabilize the calculation. Satisfactory results were 
obtained, as shown in Figures 5 and 6, when turn- 
ing off the artificial viscosity only in the near-wall 
area while having DIS4=0.6 elsewhere. Regarding 
the computational cost, 1000 iterations on the Cray 
YMP computer took 177 seconds for the wall func- 
tion approach and 262 seconds for the low Reynolds 
number approach. 

5. Conclusions 

The compressible wall functions based on the 
Van Driest transformation have been implemented 
in the turbulence module developed for the NPARC 
code. The details of the numerical implementation 
are reported. 

The module is validated against the incompress- 
ible flow over a flat plate. Very good results have 
been obtained. It is found that the unrealistic artifi- 
cial viscosity near the wall has a big influence on the 
solution of the wall function approach. The simplest 
way to eliminate this influence, as used in this work, 
is to turn off the artificial viscosity in the near-wall 
region. The module has to undergo further tests be- 
fore fully establishing its computational capability. 

This work is only restricted to the simple form 
of the wall functions. Refined forms have been avail- 
able which account for pressure gradient effects and 


mass flow through the wall or the effects of the vis- 
cous sublayer. Their implementation will be a sub- 
ject of further development of the module. 
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Fig. 1 A typical control volume centered at node C 
and its neighbouring nodes 



Fig.2 Geometry of solution field 
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Fig. 5 U- velocity profile 
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Fig. 4 Influence of artificial viscosity 
on wall friction coefficient 
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Fig. 6 Wall friction coefficient 
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